//#pragma OPENCL EXTENSION all : enable

#define ALIGN     16
#define ROUND(n)  (((size_t)(n)+ALIGN-1)&(~(size_t)(ALIGN-1)))

#define SPEED_MIN (float) -0.2
#define SPEED_MAX (float) 0.2

#define SPEED_ASSIGN(x, y) atomic_xchg(&speed[x], clamp(y, SPEED_MIN, SPEED_MAX));
#define COPY_VECTOR(dest, vect, i) { dest.x = vect[i]; 	dest.y = vect[i + ROUND(N)]; dest.z = vect[i + 2*ROUND(N)]; }

void add_speed(__global float *speed, int index, float dspeed)
{
	SPEED_ASSIGN(index, speed[index] + dspeed)
}

// This kernel is executed with a very high number of threads = total number of vertices * 3
__kernel
void eating(__global float *vbo, float dy, unsigned vertices_per_atom)
{
   int index = get_global_id(0);
   int tpa = 3 * vertices_per_atom; // threads per atom
   // ghost or pacman?
   float enable = (index / tpa) % 2 ? 0.0f : 1.0f;
   // y coordinate?
   float is_y = ((index % 3) == 1) ? 1.0f : 0.0f;
   // vertex belongs to upper part of sphere?
   float is_up = (index % tpa >= tpa/2) ? 1.0f : -1.0f;

   vbo[index] += dy * is_y * is_up * enable;
}

__kernel
void move_vertices(__global float *vbo, __global float *speed, unsigned vertices_per_atom)
{
	int index = get_global_id(0);
	int tpa = 3 * vertices_per_atom; 
	
	int x = index/tpa;
	int offset = ROUND(get_global_size(0)/tpa);
	
	vbo[index] += speed[(index%3)*offset+x];
}

__kernel
void update_position(__global float *pos, __global float *speed, unsigned N)
{

	int id = get_global_id(0);
	int index;

	index = id%N + id/N*ROUND(N);

	if(id % ROUND(N)<N)
	{
		pos[id] += speed[id];
	}
}

__kernel
void border_collision(__global float *pos, __global float *speed, __constant float *min, __constant float *max, float radius, unsigned N)
{
	int id = get_global_id(0);
	int index  =  id / ROUND(N);	
	
	if((pos[id] - radius < min[index] && speed[id] < 0) 
	|| (pos[id] + radius > max[index] && speed[id] > 0))
	{
		speed[id] *= -1.0;
	}
}

void collide_compute(__global float *pos, __global float *speed, float radius, float3 m1,  int id, int N, int p)
{
	float3 m2;
	COPY_VECTOR(m2, pos, p)
	float dist = distance(m1, m2);
	
	if(dist < 2*radius) //p-ê faire un || sur chaque axe ?
	{
		float3 v1, v2, M1, M2, M3;
		M1 = normalize(m2 - m1); 

		//précalculer 1/(1+M1.x)
		float delta =  1 / (1 + M1.x);
		M2.x = -M1.y;
		M2.y = M1.x + M1.z*M1.z * delta;
		M2.z = -M1.y * M1.z * delta;
		
		M3.x = -M1.z;
		M3.y = M2.z;
		M3.z = M1.x + M1.y*M1.y * delta;
		
		COPY_VECTOR(v1, speed, id)
		COPY_VECTOR(v2, speed, p)
		
		float3 Var = {dot(M1, v2), dot(M2, v1), dot(M3, v1)};
		float3 Vbr = {dot(M1, v1), dot(M2, v2), dot(M3, v2)};
	
		SPEED_ASSIGN(id,                M1.x * Var.x + M2.x * Var.y + M3.x * Var.z)
		SPEED_ASSIGN(id + ROUND(N),     M1.y * Var.x + M2.y * Var.y + M3.y * Var.z)
		SPEED_ASSIGN(id + 2 * ROUND(N), M1.z * Var.x + M2.z * Var.y + M3.z * Var.z)
		
		SPEED_ASSIGN(p,                 M1.x * Vbr.x + M2.x * Vbr.y + M3.x * Vbr.z)
		SPEED_ASSIGN(p + ROUND(N),      M1.y * Vbr.x + M2.y * Vbr.y + M3.y * Vbr.z)
		SPEED_ASSIGN(p + 2 * ROUND(N),  M1.z * Vbr.x + M2.z * Vbr.y + M3.z * Vbr.z)
					
		if(dist < radius 
		&& speed[p] < dist
		&& speed[p + ROUND(N)] < dist 
		&& speed[p + 2 * ROUND(N)] < dist)
		{	
			SPEED_ASSIGN(p, ((speed[p] < 0) ? -1 : 1) * dist) 
		} 
	}
}

__kernel
void atom_collision(__global float *pos, __global float *speed, float radius, int N)
{
	int id = get_global_id(0);
	int NumThreads = get_global_size(0);
	float3 m1;
	
	
	/* Methode 1 */
	/*
	COPY_VECTOR(m1, pos, id)
	for(int p = id+1; p<N; p++)
	{
		collide_compute(pos, speed, radius, m1, id, N, p);
	}
	*/
	
	/* Méthode 2 */
	COPY_VECTOR(m1, pos, id)
	for(int p = id+1; p < N ; ++p)
	{
		collide_compute(pos, speed, radius, m1, id, N, p);
		
		if(p >= N - id && id != NumThreads - 1)
		{
			collide_compute(pos, speed, radius, m1, N - id - 1 , N, p);
		}
		
		//~ if(N%2 && id == NumThreads - 1 && p > id + 1) // cas impair
		//~ {
			//~ collide_compute(pos, speed, radius, m1, id+1, N, p);
		//~ }
	}
}

__kernel
void gravity(__global float *pos, __global float *speed, float g)
{
	int id = get_global_id(0);
	int N = get_global_size(0);
	
	int	ycoord = id + N;

	if(speed[ycoord] >= 0) // on monte
	{
		add_speed(speed, ycoord,  -2 * g); //un peu d'inertie
	}
	else
	{
		add_speed(speed, ycoord, -g);
	}
}

__kernel
void lennard_jones(__global float *pos, __global float *speed, float r0)
{
	int id = get_global_id(0);
	int N = get_global_size(0);
	int i = 0;
	float3 V, V2, m1, m2;
	float r, Ep;
	// pour chaque atome, on calcule le vecteur normé envers tous les autres atomes

	COPY_VECTOR(m1, pos, id)

	for(i = id+1; i < N; ++i)
	{
		COPY_VECTOR(m2, pos, i)
	
		V = normalize(m1 - m2);
		r = distance(m1, m2);
		
		// s = 2 r0
		// formule : 24 * eps / sigma * ((2s/r)^13 - (s/r)^7)
		V *= 24 / r0 * (pown(3.4  * r0 / r, 13) - pown(2 * r0 / r, 7));

		add_speed(speed, id, V.x);
		add_speed(speed, id + ROUND(N), V.y);
		add_speed(speed, id + 2 * ROUND(N), V.z);

		add_speed(speed, i, -V.x);
		add_speed(speed, i + ROUND(N), -V.y);
		add_speed(speed, i + 2 * ROUND(N), -V.z);
	}
}
